Longitudinal and survival data modelisation

Initialisation et génération des données

## [1] "number of cores available = 1"

#Phi[1] ; eta = valeur de fin  Phi[2] = valeur du noeud  Phi[3] = echelle
m <- function(t, eta, phi) (phi[,1] + eta)/(1+exp((phi[,2]-t)/phi[,3]))
#=======================================#
param <- list(sigma2 = 0.05,
              rho2 = 0.1,
              mu = c(5,90,5),
              omega2 = c(0.5,0.1,0.01),
              #Survival data,
              nu2 = 0.5,
              a = 90,
              b = 50,
              alpha = 7,
              beta = 10)

#=======================================#
t <- seq(60,120, length.out = 10) #value of times

set.seed(18031997)
data <- JLS_data(G = 8, ng = 4, time = t, fct = m, param = param, linkfct = function(t,eta,phi) phi[,2]/phi[,3])

getDim(data)

##   G  ng   N   n  F. 
##   8   4  32 320   3

Y <- data$obs

Statistique exhaustive

\[\log f(Y, Z ; \theta) = \langle \Phi(\theta) ; S(Y,Z) \rangle - \psi(\theta)\]

sigma2_a <- sigma2_b <- sigma2_alpha <- 0.1

Phi <- fct_vector(function(sigma2, rho2, mu, omega2, nu2, a, b, alpha, beta) {
  c(- n/(2*sigma2), -N/(2*rho2),               #1, 2
        - G/(2*omega2), mu,                    #3, 4
        - G/(2*nu2),                               #5
          a/sigma2_a,         -1/(2*sigma2_a),     #6, 7
          b/sigma2_b        , -1/(2*sigma2_b),     #8, 9
          alpha/sigma2_alpha, -1/(2*sigma2_alpha), #10, 11
          beta ) },#12
  dim = c(1,1,F., F., rep(1,8)) )$eval

zeta <- function(beta, eta, phi, gamma, a,b, alpha)
{
  lbd <- function(t,g) t^{b-1} * exp(alpha*data@linkfct(t, eta[g], matrix(phi[g,], nrow = 1)))
  P1 <- 1:nrow(data@survival) %>% sapply(function(i) lbd(data@survival$obs[i], data@survival$gen[i]) )
    # sapply(function(i) integrate(function(t) lbd(t,data@survival$gen[i]), 0, data@survival$obs[i])$value )

  P2 <- exp(beta*data@survival$U + as(rep(gamma, each = ng), 'matrix'))
  
  beta*sum(data@survival$U) -  b*a^-b * sum(P2*P1)
}
zeta.der <- function(beta, eta, phi, gamma, a, b, alpha)
{
  lbd <- function(t,g) t^{b-1} * exp(alpha*data@linkfct(t, eta[g], matrix(phi[g,], nrow = 1)))
  P1 <- 1:nrow(data@survival) %>% sapply(function(i) data@survival$obs[i]/b*lbd(data@survival$obs[i], data@survival$gen[i]) )
    # sapply(function(i) integrate(function(t) lbd(t,data@survival$gen[i]), 0, data@survival$obs[i])$value )

  P2 <- data@survival$U*exp(beta*data@survival$U + as(rep(gamma, each = ng), 'matrix'))
  
  sum(data@survival$U) -  b*a^-b * sum(P2*P1)
}

S <- fct_vector(function(eta, phi, ...) mean((Y - get_obs(data, eta = eta, phi = phi) )^2 ),
                function(eta, ...)      mean(eta^2),           #2
                function(phi, ...)      apply(phi^2, 2, mean), #3
                function(phi, ...)      apply(phi, 2, mean),   #4
                function(gamma, ...)    mean(gamma^2),         #5
                function(a, ...)        a,                     #6
                function(a, ...)        a^2,                   #7
                function(b, ...)        b,                     #8
                function(b, ...)        b^2,                   #9
                function(alpha, ...)    alpha,                 #10
                function(alpha, ...)    alpha^2,               #11
                function(eta, phi, gamma, a, b, alpha)                        #12
                  uniroot(function(beta)zeta.der(beta, eta, phi, gamma, a,b,alpha), 
                          lower = -1000, upper = 1000)$root,#, tol = 1e-3)$root,
                dim = c(1,1,F., F., rep(1,8)) ) 

h <- fct_vector(function(a,b, ...) N*log(b*a^-b),
                function(b, ...) (b-1)*sum(log(t)),
                function(gamma,...) ng*sum(gamma),
                function(eta, phi, ...) sum(data@fct(1, eta, apply(phi, 2, rep, ng) )))

S[12](data@eta, data@phi, data@gamma, param$a, param$b, param$alpha)

## [1] 4.708159
## attr(,"12")
## [1] 1

Metropolis Hastings

Vraisemblance

loglik.phi <- function(x, eta, gamma, a,b,alpha,Phi){
  id <- c(1,3,4)
  xp <- setoffset(x, Phi%a%4)
  sum(Phi%a%1*S$eval(eta = eta, phi = xp, i = 1) + Phi%a%3 * S$eval(phi = x, i = 3) ) +
    zeta(Phi%a%12, eta, xp, gamma, a, b, alpha)
}
loglik.eta <- function(x, phi, gamma, a,b,alpha, Phi){ 
  id <- c(1,2)
  sum( Phi%a%id * S$eval(eta = x, phi = phi, i = id) ) +
    zeta(Phi%a%12, x, phi, gamma, a, b, alpha)
}

loglik.gamma <- function(x, eta, phi, a, b, alpha, Phi){
  Phi%a%5  * sum(x^2) + ng*sum(x) + h$eval(gamma = x, i = 3) +
    zeta(Phi%a%12, eta, phi, x, a, b, alpha) #Correction fait à la va vite
}

#c'est faux à refaire !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
loglik.a <- function(x, b, Phi) sum( Phi%a%6:7 * S$eval(a = x, b = b, i = 6:7) ) +  h$eval(a = x, b = b, i = 1)
loglik.b <- function(x, a, Phi) sum( Phi%a%8:9 * S$eval(b = x, i = 8:9) ) +  sum(h$eval(a = a, b = x, i = c(1,2)))
loglik.alpha <-  function(x, eta, phi, gamma, a,b,Phi){
  id <- 10:11
  sum( Phi%a%id * S$eval(alpha = x, i = id) ) +
    zeta(Phi%a%12, eta, phi, gamma, a, b, x) +
    x*h$eval(eta = eta, phi = phi, i = 4)
}

SAEM

initialisation

# ---  Initialisation des paramètres --- #
para <- param %>% sapply(function(x) x* runif(1, 1.1,1.4))
# para$rho2 = 0.2 ; para$omega2 <- rep(.1,3)

# --- Initialisation des chaines MC : Z_0 --- #
Z <- getLatente(data, format = 'list') #list()
Z$eta = list( matrix(rep(0, G*ng), ncol = 1) )
Z$phi = list(  matrix(rep(para$mu, G), nrow = F.) %>% t  )

Z$gamma = list(matrix(rep(0, G), ncol = 1))

Z$a = list(matrix(para$a))
Z$b = list(matrix(para$b))
Z$alpha = list(matrix(para$alpha))

Etape simulation et maximisation du SEAM

sim <- function(niter, h, Phih, eta, phi, gamma, a, b, alpha, verbatim = F)
{
  eta <- list(MH_High_Dim_future(niter, eta[[1]], sd = sd.eta, loglik.eta,
                                 phi[[1]], gamma[[1]], a[[1]], b[[1]], alpha[[1]],
                                 Phih, cores = ncores, verbatim = verbatim ))

  phi <- list(MH_Gibbs_Sampler_future(niter, setoffset(phi[[1]],-Phih%a%4), sd.phi, loglik.phi,
                                      eta[[1]], gamma[[1]], a[[1]], b[[1]], alpha[[1]],
                                      Phih, cores = 1, verbatim = verbatim )) %>%
    lapply(setoffset, Phih%a%4)
  
  gamma <- list(MH_High_Dim_future(niter, gamma[[1]], sd = .03, loglik.gamma,
                                   eta[[1]], phi[[1]], a[[1]], b[[1]], alpha[[1]],
                                   Phih, cores = 1, verbatim = verbatim ))

  a <- list(MH_High_Dim_future(niter, a[[1]],     sd = .5, loglik.a, b[[1]],
                               Phih, cores = 1, verbatim = verbatim ))
  b <- list(MH_High_Dim_future(niter, b[[1]],     sd = .5, loglik.b, a[[1]],
                               Phih, cores = 1, verbatim = verbatim ))
  
  alpha <- list( MH_High_Dim_future(niter, alpha[[1]], sd = 1, loglik.alpha,
                                   eta[[1]], phi[[1]], gamma[[1]], a[[1]], b[[1]],
                                    Phih, cores = 1, verbatim = verbatim ))

  list(eta = eta, phi = phi, gamma = gamma, a = a, b = b, alpha = alpha)
}

#res <- simulation_test(sim, Phi, param, 1000, Z) %>% lapply(function(z)plot(z[[1]], nrow = 2))
#names(res) %>% sapply(function(var) assign(var, res[[var]][[1]],envir = .GlobalEnv))
# oracle(maxi, S, data, a = param$a, b = param$b, alpha = param$alpha)

maxi <- function(S)
{
  list(sigma2 = param$sigma2,#S%a%1,
       rho2 =   param$rho2,#S%a%2,
       mu =     param$mu,#S%a%4,
       omega2 = param$omega2,#S%a%3 - (S%a%4)^2,
       nu2 =    S%a%5,
       a =      S%a%6,  b =    S%a%8,
       alpha =  S%a%10, beta = param$beta)#S%a%12 )
}

Resultats

niter <- 20
MH.iter <- 10

sd.eta <- .3
sd.phi <- c(.05, .3, .05)

gg <- simulation_test(sim, Phi, param, 250, Z) %>% lapply(function(z)plot(z[[1]], nrow = 2))

res <- SAEM(niter, MH.iter, para, Phi, S$eval, Z, sim, maxi, eps = 1e-2,
            burnin = 150, coef.burnin = 2/3, verbatim = 2)
saveRDS(res, 'saem.rds')

## [1] "SAEM execution time = 01min 59sec"

## $plot_parameter

## 
## $plot_MCMC

## 
## $plot_acceptation
## TableGrob (3 x 2) "arrange": 6 grobs
##   z     cells    name           grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
## 3 3 (2-2,1-1) arrange gtable[layout]
## 4 4 (2-2,2-2) arrange gtable[layout]
## 5 5 (3-3,1-1) arrange gtable[layout]
## 6 6 (3-3,2-2) arrange gtable[layout]

[1] “SAEM execution time = 01min 59sec”

Result of the SAEM-MCMC

	sigma2	rho2	mu.1	mu.2	mu.3	omega2.1	omega2.2	omega2.3	nu2	a	b	alpha	beta
Real value	0.05	0.1	5	90	5	0.5	0.1	0.01	0.2899	90.0000	50.0000	7.0000	10
Estimated value	0.05	0.1	5	90	5	0.5	0.1	0.01	0.0000	74.0437	30.0973	2.9938	10
Rrmse	0.00	0.0	0	0	0	0.0	0.0	0.00	1.0000	0.1773	0.3981	0.5723	0